home *** CD-ROM | disk | FTP | other *** search
/ TPUG - Toronto PET Users Group / TPUG Users Group CD / TPUG Users Group CD.iso / AMIGA / AMICUS / AMICUS18.ADF / Progs / StarProbe / SPEVAL.C < prev    next >
C/C++ Source or Header  |  1989-01-27  |  12KB  |  495 lines
  1. /****************************************************************************
  2.  ***                                                                      ***
  3.  ***      S T A R   P R O B E        E V A L U A T E                      ***
  4.  ***                                                                      ***
  5.  ****************************************************************************/
  6. /*****
  7.  ***  This module evaluates all physical parameters at the current mass
  8.  ***  position.  All physics/calculus emanates from here.
  9.  *****/
  10.  
  11. #include "stdio.h"
  12. #include "math.h"
  13. #include "limits.h"
  14. #include "spmac.h"
  15. #include "spref.h"
  16.  
  17. int ixstop,ixstart;
  18. double m2,init_p,init_t,init_l,init_r,massaccum;
  19.  
  20. /*****
  21.  ***   M e a n   M o l e c u l a r   W e i g h t
  22.  *****/
  23.  
  24. double meanmolweight(h_frac,he_frac,m_frac)
  25. double h_frac,    /* fraction of hydrogen */
  26.        he_frac,   /* fraction of helium */
  27.        m_frac;    /* fraction of metals */
  28. {                 /* h+he+m frac must = 1.0 */
  29. iam(80);
  30. static double val;
  31.  
  32.  block(in,"meanmolweight",(h_frac+he_frac+m_frac));
  33.  assert0((h_frac+he_frac+m_frac),==,1.0);
  34.  
  35.  val = (1.0/((2*h_frac)+(0.75*he_frac)+(0.5*m_frac)));
  36.  
  37.  block(out,"meanmolweight",val);
  38.  return(val);
  39. }
  40.  
  41.  
  42. /*****
  43.  ***   E v a l D e n s i t y
  44.  ***
  45.  ***   Ref: Motz 12.11, pg. 282
  46.  *****/
  47.  
  48. double eval_density(p,t)
  49. double p,t;
  50. {
  51. iam(81);
  52. static double val;
  53. static double mpk = MP/K;
  54.  block(in,"eval_density",p);
  55.  
  56.  assert1(t,>,0.0);
  57.  
  58.  val = u*mpk*p / t;
  59.  
  60.  assert1(val,<,500.0);
  61.  
  62.  block(out,"eval_density",val);
  63.  return(val);
  64. }
  65.  
  66. /*****
  67.  ***   E v a l E n e r g y R a t e
  68.  ***
  69.  ***   Ref: Motz pp. 344-346 (13.58)
  70.  *****/
  71. double eval_energy_rate(d,t)
  72. double d,t;
  73. {
  74. iam(41);
  75. double val,t6,e_pp,e_cn,p,it;
  76. int pp_ix, cn_ix;
  77.  block(in,"eval-energy",t);
  78.  
  79.  it = pow((1.0e6/t),0.33333);
  80.  
  81.  e_pp = 2.36e6 *d*X*X*it*it*exp(-33.8*it);
  82.  
  83.  e_cn = 0.786e28 *d*X*Z*it*it*exp(-152.3*it)*exp(7.25/t);
  84.  
  85.  val = e_pp + e_cn;
  86.  
  87.  block(mid,"eval-energy e_pp",e_pp);
  88.  block(mid,"eval-energy e_cn",e_cn);
  89.  block(out,"eval-energy",val);
  90.  return(val);
  91. }
  92.  
  93. /*****
  94.  ***   E v a l O p a c i t y
  95.  ***
  96.  ***   Ref: Motz pp. 165-166
  97.  *****/
  98.  
  99. double eval_opacity(p,t)
  100. double p,t;
  101. {
  102.  iam(31);
  103.  static double val,e1,e2,d1,d2,t2,t3,t4,d,t6,d3,d4,d5;
  104.  block(in,"eval-opacity",0.0);
  105.  
  106.  t6 = t * 1.0e-6;
  107.  d  = eval_density(p,t);
  108.  t2 = t6*t6; t3 = t2*t6; t4 = t3*t6;
  109.  d1 = 6.5e4 * Z * (1.0-(0.05*t6))/(t2+(2.5*t4));
  110.  d1*= exp(-7.75*d*(1.0+X)/t3);
  111.  d2 = -abs((0.667-2.873*t6));
  112.  d2 = 1.0 + 5.5*exp(d2);
  113.  d3 = X/((250.0*t4)-t2);
  114.  d4 = Y/(250.0*t4);
  115.  d5 = 4.15e4 * (d3 + (d4*d2));
  116.  e1 = d1 + (d5 * exp(-2.0*d*(1.0+X)/t3));
  117.  
  118.  e2 = (225.0-(115.0*X)-(190.0*Y))/pow(t6,3.5);
  119.  
  120.  if (t6<10.0){
  121.    val = (0.19*(1.0+X))+(e1*d*(1.0+X));
  122.    val *= 0.715;
  123.    }
  124.  else if (t6<=20.0){
  125.    val = (0.19*(1.0+X))+(e1*d*(1.0+X))+
  126.          (d*(1.0-0.1*t6)*((1.0+X)*e1-e2));
  127.    val *= 0.9;
  128.    }
  129.  else
  130.    val = (0.19*(1.0+X))+e2*d;
  131.  
  132.  assert1(val,>,0.0);
  133.  block(out,"eval-opacity",val);
  134.  return(val);
  135. }
  136.  
  137. /*****
  138.  ***   R a d i a t i v e
  139.  ***
  140.  ***  Ref: Motz pg. 311
  141.  *****/
  142.  
  143. int radiative(m,p,t,l,r)
  144. double m,p,t,l,r;
  145. {
  146.  iam(82);
  147.  int val;
  148.  block(in,"radiative",t);
  149.  
  150.  assert2(p,>,0.0);
  151.  assert2(l,>,0.0);
  152.  assert2(m,>,0.0);
  153.  if (pow(t,5.0)/(eval_opacity(p,t)*p*l/m) <= 1.0e10)
  154.    val = 0;
  155.  else
  156.    val = 1;
  157.  
  158.  block(out,"radiative",(double)val);
  159.  return(val);
  160. }
  161.  
  162. /*****
  163.  ***   P r e s s u r e   F u n c t i o n
  164.  ***
  165.  ***   Ref: Motz pg. 292 or Clayton Eq 6-14
  166.  *****/
  167.  
  168. double press_func(m,p,t,l,r)
  169. double m,p,t,l,r;
  170. {
  171. iam(20);
  172. static double val;
  173.  block(in,"press-func",m);
  174.  
  175.  assert2(r,>,0.0);
  176.  
  177.  val=((-G*m)/(4.0*PI*r*r*r*r));
  178.  
  179.  block(out,"press-func",val);
  180.  return(val);
  181. }
  182.  
  183. /*****
  184.  ***   T e m p e r a t u r e   F u n c t i o n
  185.  ***
  186.  ***   Ref: Clayton pg. 443, Motz pp. 252-253
  187.  *****/
  188.  
  189. double temp_func(m,p,t,l,r)
  190. double m,p,t,l,r;
  191. {
  192. iam(30);
  193. double tau2, val, opcty;
  194.  block(in,"temp-func",t);
  195.  
  196.  if (radiative(m,p,t,l,r)){
  197.    opcty = eval_opacity(p,t);
  198.    assert2(t,>,0.0);
  199.    assert2(r,>,10000.0);
  200.    val = ((-3.0*opcty*l)/
  201.           (4.0*A*C*t*t*t*16.0*PI_SQ*r*r*r*r));
  202.    }
  203.  else{
  204.    tau2 = 1.0 + ((4.0-3.0*B)*(poly-1.0))/(B*B+3.0*(poly-1.0)*(1.0-B)*(4.0+B));
  205.    block(mid,"temp-func tau2",tau2);
  206.    val = ((tau2-1.0)*t/(tau2*p)) * press_func(m,p,t,l,r);;
  207.    }
  208.  
  209.  block(out,"temp-func",val);
  210.  return(val);
  211. }
  212.  
  213. /*****
  214.  ***   L u m i n o s i t y   F u n c t i o n
  215.  ***
  216.  ***   Ref: Motz pg. 292
  217.  *****/
  218.  
  219. double lum_func(m,p,t,l,r)
  220. double m,p,t,l,r;
  221. {
  222. iam(40);
  223. static double val,d;
  224.  block(in,"lum-func",l);
  225.  
  226.  d = eval_density(p,t);
  227.  val = eval_energy_rate(d,t);
  228.  assert2(val,>=,0.0);
  229.  
  230.  block(out,"lum-func",val);
  231.  return(val);
  232. }
  233.  
  234. /*****
  235.  ***   R a d i u s   F u n c t i o n
  236.  ***
  237.  ***   Ref: Clayton Eq 6-13
  238.  *****/
  239.  
  240. double radius_func(m,p,t,l,r)
  241. double m,p,t,l,r;
  242. {
  243. iam(50);
  244. static double val;
  245.  block(in,"radius-func",r);
  246.  
  247.  assert2(r,>,0.0);
  248.  val = 1.0/(4.0*PI*r*r*eval_density(p,t));
  249.  
  250.  block(out,"radius-func",val);
  251.  return(val);
  252. }
  253.  
  254.  
  255. /*****
  256.  ***   F u n c
  257.  ***
  258.  ***   Calls the correct differential function for the integrator.
  259.  *****/
  260.  
  261. double func(i,m,p,t,l,r)
  262. int i;
  263. double m,p,t,l,r;
  264. {
  265. iam(13);
  266. static double val;
  267.  block(in,"func",(double)i);
  268.  block(mid,"func m",m);
  269.  block(mid,"func p",p);
  270.  block(mid,"func t",t);
  271.  block(mid,"func l",l);
  272.  block(mid,"func r",r);
  273.  
  274.  switch (i){
  275.    case 0: val = press_func(m,p,t,l,r);
  276.            break;
  277.    case 1: val = temp_func(m,p,t,l,r);
  278.            break;
  279.    case 2: val = lum_func(m,p,t,l,r);
  280.            break;
  281.    case 3: val = radius_func(m,p,t,l,r);
  282.            break;
  283.    default:assert_error(proc_num,(double)i,4.0,0.0,0.0,0.0,0.0,0.0);
  284.    }
  285.  
  286.  block(out,"func",val);
  287.  return (val);
  288. }
  289.  
  290. /*****
  291.  ***   I n t e g r a t e
  292.  ***
  293.  *** Solves for (P,T,L,R) using Runge_Kotta algorithm. 
  294.  ***
  295.  *** Ref: Clayton pp. 448-450
  296.  *****/
  297.  
  298. void integrate(ixstrt,ixstop,ixinc,p,t,l,r)
  299. int ixstrt,ixstop,ixinc;/* array index control */
  300. double p,t,l,r;
  301. {
  302.  iam(12);
  303.  double k1[4], k2[4], k3[4], k4[4], h, h2, mh, mh2;
  304.  double delta_p,delta_t,delta_l,delta_r;
  305.  int i,f;
  306.  block(in,"integrate",0.0);
  307.  
  308.  i = ixstrt;
  309.  while (1){
  310.    block(mid,"integrate ***************************",(double) i);
  311.    assert0(i,>=,0);
  312.    assert0(i,<,gradsteps);
  313.    massinc=massgrad[i]*ixinc*mass;
  314.    h = massinc; h2 = h/2.0;
  315.    mh = m + h;  mh2 = m + h2;
  316.  
  317.    block(mid,"integrate m",m);
  318.    block(mid,"integrate h",h);
  319.    block(mid,"integrate p",p);
  320.    block(mid,"integrate t",t);
  321.    block(mid,"integrate l",l);
  322.    block(mid,"integrate r",r);
  323.  
  324.    block(mid,"integrate *** K1 ****",0.0);
  325.    for (f=0;f<4;f++)
  326.      k1[f] = func(f,m,p,t,l,r);
  327.  
  328.    block(mid,"integrate *** K2 ****",0.0);
  329.    for (f=0;f<4;f++)
  330.      k2[f] = func(f,mh2,p+h2*k1[0],t+h2*k1[1],l+h2*k1[2],r+h2*k1[3]);
  331.  
  332.    block(mid,"integrate *** K3 ****",0.0);
  333.    for (f=0;f<4;f++)
  334.      k3[f] = func(f,mh2,p+h2*k2[0],t+h2*k2[1],l+h2*k2[2],r+h2*k2[3]);
  335.  
  336.    block(mid,"integrate *** K4 ****",0.0);
  337.    for (f=0;f<4;f++)
  338.      k4[f] = func(f,mh,p+h*k3[0],t+h*k3[1],l+h*k3[2],r+h*k3[3]);
  339.         
  340.    delta_p=(h*(k1[0]+k2[0]+k2[0]+k3[0]+k3[0]+k4[0])/6.0);
  341.    delta_t=(h*(k1[1]+k2[1]+k2[1]+k3[1]+k3[1]+k4[1])/6.0);
  342.    delta_l=(h*(k1[2]+k2[2]+k2[2]+k3[2]+k3[2]+k4[2])/6.0);
  343.    delta_r=(h*(k1[3]+k2[3]+k2[3]+k3[3]+k3[3]+k4[3])/6.0);
  344.  
  345.    p += delta_p;
  346.    t += delta_t;
  347.    l += delta_l;
  348.    r += delta_r;
  349.  
  350.    pressure = p;
  351.    temp = t;
  352.    lum = l;
  353.    radius = r;
  354.    density = eval_density(p,t);
  355.    opacity = eval_opacity(p,t);
  356.    pressgrad[i] = p;
  357.    tempgrad[i] = t;
  358.    lumgrad[i] = l;
  359.    radiusgrad[i] = r;
  360.    densgrad[i] = density;
  361.    opctygrad[i] = (radiative(m,p,t,l,r)) ? opacity : 0.0;
  362.    m+=massinc;
  363.    tmassgrad[i] = m;
  364.    energygrad[i] = eval_energy_rate(density,t);
  365.  
  366.    block(mid,"integrate delta_p=",delta_p);
  367.    block(mid,"integrate delta_t=",delta_t);
  368.    block(mid,"integrate delta_l=",delta_l);
  369.    block(mid,"integrate delta_r=",delta_r);
  370.  
  371.    if (debug > 0) report_detail();
  372.    if (i==ixstop) break;
  373.    i += ixinc;
  374.    }
  375.  
  376.  results[0] [rix] = p;
  377.  results[1] [rix] = t;
  378.  results[2] [rix] = l;
  379.  results[3] [rix] = r;
  380.  
  381.  block(mid,"integrate ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! !",0.0);
  382.  block(out,"integrate",0.0);
  383. }
  384.  
  385. /*****
  386.  ***   C a l l   I n t e g r a t o r
  387.  ***
  388.  ***   Sets up the global variables and parameters for the Integrate func.
  389.  *****/
  390.  
  391. void call_integrator(dir,ix,p_mult,t_mult,l_mult,r_mult)
  392. int dir,ix;
  393. double p_mult,t_mult,l_mult,r_mult;
  394. {
  395. iam(14);
  396.  block(in,"call-integrator",dir);
  397.  block(mid,"call-integrator",ix);
  398.  
  399.  if (dir==in) { /* inward from surface to fitting point */
  400.    ixstart = gradsteps - 1;
  401.    /*init_p = 0.667*G*SM/(SR*SR*radratio*radratio*10.0);*/
  402.    init_p = effpress * p_mult;
  403.    init_t = efftemp * t_mult;
  404.    init_r = SR*radratio * r_mult;
  405.    /*init_l = 4.0*PI*init_r*init_r*SB*pow(efftemp,4.0);*/
  406.    init_l = SL*lumratio * l_mult;
  407.  
  408.    rix=ix;m=mass;
  409.    pressure=init_p;temp=init_t;lum=init_l;radius=init_r;
  410.    density=eval_density(init_p,init_t);opacity=eval_opacity(init_p,init_t);
  411.    if (trace_me) report_detail();
  412.    integrate(ixstart,ixstop+1,-1,init_p,init_t,init_l,init_r);
  413.    }
  414.  else { /*outward from center to fitting point */
  415.    init_p = corepress * p_mult;
  416.    init_t = coretemp * t_mult;
  417.    init_l = l_mult * SL * lumratio / 1.0e8;
  418.    init_r = r_mult * pow((3.0*mass*massgrad[0]/(4.0*PI*coredensity)),0.33333);
  419.  
  420.    rix=ix;m=mass*massgrad[0];
  421.    pressure=init_p;temp=init_t;lum=init_l;radius=init_r;
  422.    density=eval_density(init_p,init_t);opacity=eval_opacity(init_p,init_t);
  423.    pressgrad[0] = init_p;
  424.    tempgrad[0] = init_t;
  425.    lumgrad[0] = init_l;
  426.    radiusgrad[0] = init_r;
  427.    densgrad[0] = density;
  428.    opctygrad[0] = (radiative(m,init_p,init_t,init_l,init_r)) ? opacity : 0.0;
  429.    tmassgrad[0] = m;
  430.    energygrad[0] = eval_energy_rate(density,init_t);
  431.    if (trace_me) report_detail();
  432.    integrate(1,ixstop,1,init_p,init_t,init_l,init_r);
  433.    }
  434.  block(out,"call-integrator",0.0);
  435. }
  436.  
  437. /*****
  438.  ***   E v a l u a t e
  439.  ***
  440.  ***   Calls Integrate 6 times for each trial integration: 
  441.  *** 3 from the center, 3 from the surface; each ending when
  442.  *** m=mass/2. The results are then massaged to get new starting
  443.  *** values for the next trial. The trials stop when the results
  444.  *** become consistent.
  445.  ***
  446.  ***   The initial values come from Set_Boundary_Values.
  447.  *****/
  448.  
  449. void evaluate()
  450. {
  451.  iam(10);
  452.  int tries;
  453.  block(in,"evaluate",0.0);
  454.  
  455.  massaccum=0.0;
  456.  m2=mass*fitptmass;
  457.  ixstop=-1;
  458.  do {         /* find mass midpoint */
  459.    ixstop++;
  460.    massaccum+=massgrad[ixstop]*mass;
  461.    } while (massaccum<m2);
  462.  assert0(massaccum,>,0.0);
  463.  assert0(massaccum,<=,mass); 
  464.  block(mid,"evaluate ixstop=",(double) ixstop);
  465.  
  466.  tries = 0;
  467.  for (;;) {
  468.    tries++;
  469.    report_global();
  470.    printf("Starting trial integration #%d\n",tries);
  471.  
  472.    call_integrator(out,0,1.0,  1.0,  1.0,  1.0  );
  473.    call_integrator(in, 1,1.0,  1.0,  1.0,  1.0  );
  474.    dump_data(tries);
  475.  
  476.    if (!discontinuity()) break;
  477.    if (tries >= max_tries) break;
  478.  
  479.    if (!closer_fit(ixstop)){
  480.      printf("... integrating with variances ...\n");
  481.      call_integrator(out,2,adj_p,1.0,  1.0,  1.0  );
  482.      call_integrator(out,3,1.0,  adj_t,1.0,  1.0  );
  483.      call_integrator(in, 4,1.0,  1.0,  adj_l,1.0  );
  484.      call_integrator(in, 5,1.0,  1.0,  1.0,  adj_r);
  485.      }
  486.  
  487.    report_results();
  488.  
  489.    new_boundary_conditions();
  490.    }
  491.  
  492.  block(out,"evaluate",0.0);
  493. }
  494.  
  495.